SKELLAM

Overview

The SKELLAM function computes values for the Skellam distribution, a discrete probability distribution that describes the difference between two independent Poisson-distributed random variables. Named after the British statistician John Gordon Skellam, who first described it in 1946, this distribution is useful for modeling the difference in counts from two independent Poisson processes.

The Skellam distribution arises naturally when comparing two count-based phenomena. Given two independent Poisson random variables N_1 and N_2 with expected values \mu_1 and \mu_2 respectively, the difference K = N_1 - N_2 follows a Skellam distribution. This makes it applicable to scenarios such as analyzing the point spread in sports like soccer, hockey, or baseball where both teams score discrete points, as well as comparing photon counts in imaging applications.

The probability mass function (PMF) is given by:

p(k; \mu_1, \mu_2) = e^{-(\mu_1 + \mu_2)} \left(\frac{\mu_1}{\mu_2}\right)^{k/2} I_{|k|}(2\sqrt{\mu_1 \mu_2})

where I_k(z) is the modified Bessel function of the first kind. Unlike many discrete distributions that only take non-negative values, the Skellam distribution has support on all integers (positive, negative, and zero), reflecting the possibility that either Poisson variable could be larger.

The distribution has mean \mu_1 - \mu_2 and variance \mu_1 + \mu_2. When \mu_1 = \mu_2, the distribution is symmetric around zero. Both parameters \mu_1 and \mu_2 must be strictly positive.

This implementation uses SciPy’s scipy.stats.skellam module, which provides methods for the PMF, CDF, survival function, inverse CDF (percent point function), and statistical measures including mean, variance, and standard deviation. The optional loc parameter shifts the distribution along the integer axis.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=SKELLAM(k, mu_one, mu_two, skellam_mode, loc)

k (float, required): Value at which to evaluate the distribution (or probability for ICDF/ISF modes). Can be a scalar or 2D list.
mu_one (float, required): Mean of the first Poisson variable (must be greater than 0).
mu_two (float, required): Mean of the second Poisson variable (must be greater than 0).
skellam_mode (str, optional, default: “pmf”): Output type for the distribution calculation.
loc (float, optional, default: 0): Location parameter that shifts the distribution.

Returns (float): Distribution result (float), or error message string.

Examples

Example 1: PMF at k=2

Inputs:

k	mu_one	mu_two	skellam_mode	loc
2	5	3	pmf	0

Excel formula:

=SKELLAM(2, 5, 3, "pmf", 0)

Expected output:

0.1431

Example 2: CDF at k=2

Inputs:

k	mu_one	mu_two	skellam_mode	loc
2	5	3	cdf	0

Excel formula:

=SKELLAM(2, 5, 3, "cdf", 0)

Expected output:

0.5776

Example 3: SF at k=2

Inputs:

k	mu_one	mu_two	skellam_mode	loc
2	5	3	sf	0

Excel formula:

=SKELLAM(2, 5, 3, "sf", 0)

Expected output:

0.4224

Example 4: ICDF at probability 0.5

Inputs:

k	mu_one	mu_two	skellam_mode	loc
0.5	5	3	icdf	0

Excel formula:

=SKELLAM(0.5, 5, 3, "icdf", 0)

Expected output:

2

Example 5: Mean of distribution

Inputs:

k	mu_one	mu_two	skellam_mode	loc
0	5	3	mean	0

Excel formula:

=SKELLAM(0, 5, 3, "mean", 0)

Expected output:

2

Python Code

from scipy.stats import skellam as scipy_skellam

def skellam(k, mu_one, mu_two, skellam_mode='pmf', loc=0):
    """
    Compute Skellam distribution values using scipy.stats.skellam.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.skellam.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        k (float): Value at which to evaluate the distribution (or probability for ICDF/ISF modes). Can be a scalar or 2D list.
        mu_one (float): Mean of the first Poisson variable (must be greater than 0).
        mu_two (float): Mean of the second Poisson variable (must be greater than 0).
        skellam_mode (str, optional): Output type for the distribution calculation. Valid options: PMF, CDF, SF, ICDF, ISF, Mean, Variance, Std Dev, Median. Default is 'pmf'.
        loc (float, optional): Location parameter that shifts the distribution. Default is 0.

    Returns:
        float: Distribution result (float), or error message string.
    """
    # Validate mu_one, mu_two
    try:
        mu_one_val = float(mu_one)
        mu_two_val = float(mu_two)
    except (TypeError, ValueError):
        return "Invalid input: mu_one and mu_two must be numbers."

    if not (mu_one_val > 0 and mu_two_val > 0):
        return "Invalid input: mu_one and mu_two must be > 0."

    # Validate loc
    try:
        loc_val = float(loc)
    except (TypeError, ValueError):
        return "Invalid input: loc must be a number."
    # Validate skellam_mode
    valid_modes = {"pmf", "cdf", "sf", "icdf", "isf", "mean", "var", "std", "median"}
    if not isinstance(skellam_mode, str) or skellam_mode not in valid_modes:
        return f"Invalid input: skellam_mode must be one of {sorted(valid_modes)}."
    # Helper to process k (scalar or 2D list)
    def process_k(val):
        try:
            return float(val)
        except (TypeError, ValueError):
            return None

    # Helper to validate probability range for ICDF/ISF modes
    def validate_probability(val, mode):
        if mode in ["icdf", "isf"]:
            if val < 0 or val > 1:
                return f"Invalid input: k must be between 0 and 1 for {mode} mode."
        return None

    # Handle statistics
    if skellam_mode == "mean":
        return float(scipy_skellam.mean(mu_one_val, mu_two_val, loc=loc_val))
    if skellam_mode == "var":
        return float(scipy_skellam.var(mu_one_val, mu_two_val, loc=loc_val))
    if skellam_mode == "std":
        return float(scipy_skellam.std(mu_one_val, mu_two_val, loc=loc_val))
    if skellam_mode == "median":
        return float(scipy_skellam.median(mu_one_val, mu_two_val, loc=loc_val))

    # PMF, CDF, SF, ICDF, ISF
    def compute(val):
        kval = process_k(val)
        if kval is None:
            return "Invalid input: k must be a number."

        # Validate probability range for ICDF/ISF modes
        prob_error = validate_probability(kval, skellam_mode)
        if prob_error:
            return prob_error

        if skellam_mode == "pmf":
            return float(scipy_skellam.pmf(kval, mu_one_val, mu_two_val, loc=loc_val))
        elif skellam_mode == "cdf":
            return float(scipy_skellam.cdf(kval, mu_one_val, mu_two_val, loc=loc_val))
        elif skellam_mode == "sf":
            return float(scipy_skellam.sf(kval, mu_one_val, mu_two_val, loc=loc_val))
        elif skellam_mode == "icdf":
            return float(scipy_skellam.ppf(kval, mu_one_val, mu_two_val, loc=loc_val))
        elif skellam_mode == "isf":
            return float(scipy_skellam.isf(kval, mu_one_val, mu_two_val, loc=loc_val))
    # 2D list or scalar
    if isinstance(k, list):
        if not all(isinstance(row, list) for row in k):
            return "Invalid input: k must be a scalar or 2D list."
        result = []
        for row in k:
            result_row = []
            for val in row:
                out = compute(val)
                if isinstance(out, str):
                    return out
                result_row.append(out)
            result.append(result_row)
        return result
    else:
        return compute(k)

Overview

Excel Usage

Examples

Python Code

Online Calculator